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Abstract 

The phase behavior of a system composed of spherical particles with a 
monomodal size distribution is investigated theoretically within the context of 
the van der Waals approximation for polydisperse fluids. It is shown how the 
binodals, spinodals, cloud-point and shadow curves as well as all the (poly- 
dispersity induced) critical points can be obtained for a variety of interaction 
potentials. The polydispersity induced modifications of the phase diagram 
(even for a polydispersity index / as small as I ~ 1.01) should be observable 
in some colloidal dispersions. 



I. INTRODUCTION 



Many of the systems encountered in soft-matter physics (such as e.g. colloidal disper- 
sions!, liquid crystals!, polymeric melts!) are fluids exhibiting one or several "polydisper- 
sities" because, by the very nature of their production process, these fluids are collections 
of complex molecular objects which are not strictly identical one to another. When these 
objects can be grouped into n sets of identical objects they form an ordinary n-component 
mixture! whereas here we will focus instead on the case where the differences between these 
objects (e.g. differences in size, shape, surface charge, chemical composition, etc.) are 
distributed in an almost continuous manner in which case they can be described as "con- 
tinuous" mixtures containing infinitely many components!. In the present study we will 
consider one of the simplest situations where the system is composed of spherical molecu- 
lar objects which differ only in size. This excludes the liquid crystals and polymeric melts 
from our consideration. We will hence focus our attention on colloidal dispersions of spher- 
ical particles with a continuous size distribution^. In practice, this size distribution can 
be either monomodal (i.e. having only one maximum) or multimodal (i.e. having several 
maxima). We will again consider only the simplest case of a system with a monomodal 
size distributio Such systems can hence be considered as polydisperse generalizations of 
single component systems (whereas systems with a multimodal size distribution are general- 
izations of multi-component mixtures). When this monomodal size distribution is strongly 
peaked around its maximum one often refers to the corresponding colloidal dispersion as 
being monodisperseS. At present, such monodisperse colloidal dispersions are often used as 
testing ground for the study of various aspects of liquid state theory!. Indeed, by clever 
chemical engineering the colloidal particles can often be endowed with properties which are 
outside the range of ordinary single component systems (e.g. they can be prepared such 
as to have hard-sphere like repulsions and short-range attractions!). Nevertheless, even in 
those colloidal dispersions which are usually denoted as monodisperse there always remains 



a residual size distribution. While the influence of this (unavoidable) polydispersity on the 
behavior of a single phase is well understood!, the first-principles study of its influence on a 
phase separation process is on the contrary still under active scrutiny^. Because the theoret- 
ical study of phase separations in polydisperse systems is faced with rather heavy technical 
problems we will consider here only the simplest possible phase separation, namely a fluid- 
fluid phase transitionB. The theoretical tool used to describe this transition will again be the 
simplest one, namely the van der Waals (vdW) approximation^. It is indeed well known 
that the vdW- approximation provides us with a description of the phase separation of a 
fluid into a dilute fluid (or gaz phase) and a dense fluid (or liquid phase) which is both 
simple and physically sound. It moreover involves a (mean-field) critical point which consti- 
tutes an important feature of many phase diagrams^. The vdW-approximation is also very 
robust and can be generalized in several ways, including for polydisperse systems0'0. In 
the present study we will therefore focus our attention on the modifications brought about 
by the polydispersity to the binodal and critical point of a vdW-fluid with a monomodal 
size distribution. Such a description will of course not always be fully realistic for a given 
colloidal dispersion but we may hope that it retains the qualitative correctness which has 
characterized the widespread use of vdW-like descriptions in the past0. Since the present 
vdW-approximation is restricted to fluid phases only it will not be possible here to sort 
out which parts of the phase diagrams obtained are metastable with respect to the solid 
phases. By focusing on the fluid phases only, we are nevertheless able to show how the 
many technical problems raised by the study of phase equilibria in polydisperse systems 
can be solved, paving hereby the way for their later extensions to other phases and more 
involved descriptions. 

In section || we summarize the thermodynamic conditions of phase coexistence in a 
polydisperse system in terms of the system's free-energyBi which has here to be viewed as a 
functional of the density distribution among the species. We also deduce the thermodynamic 
stability conditions of a polydisperse system from the convexity of this free-energy functional. 
This procedure, which is intrinsic in the sense that it does not rely on taking the limit 
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of an infinitely many component system is, as far as we are aware, new. In section |I! 
we introduce the vdW free-energy in a form already used elsewherelll (for the study of 
weakly polydisperse systems) except that here we include the cases where the polydispersity 
also affects the amplitude of the interaction potential. The monomodal parent phase size 
distributions considered in this work are introduced in section [TV] while section [V] is devoted 
to the solution method to be used subsequently for obtaining the numerical solutions of the 
integral equations which, for a polydisperse system, replace the algebraic equations governing 
the phase equilibria. This method can be generalized to any system for which the free-energy 
has an excess (over ideal) part which depends on a finite number of (generalized) moments 
of the polydispersity distribution. In section [VI| we show how the thermodynamic stability 
conditions of section [II] can be reduced to closed algebraic equations whenever the excess 
free-energy is of this particular type. In section [VI 1] we present explicit phase diagrams for 
interaction potentials with and without amplitude polydispersity. The critical behavior of a 
polydisperse vdW fluid is studied further in section |V111| while the final section |D£| contains 
our conclusions. 

II. PHASE EQUILIBRIUM IN A POLYDISPERSE FLUID 

The statistical mechanical description of a polydisperse equilibrium systemlll is equiv- 
alent to a density functional theory for a system whose number density, say p(r, a) , de- 
pends besides the position variable r (assuming spherical particles) moreover on a contin- 
uous species label a (considered to be dimensionless) . Such a theory is completely de- 
termined by the knowledge of its intrinsic Helmholtz free-energy (F) per unit volume (V), 
f = F/V, which is a function of the temperature (T) and a functional of the number density, 
/ = f(T, [p]), where [p] indicates a functional dependence on p(r, er). All thermodynamic 
properties can then be deduced from f(T, [p]) by functional differentiation^!. Below we will 
consider only (uniform and isotropic) fluid phases for which p(r, a) is independent of r, hence 
p(r, a) = p{cx). 
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A. The thermodynamic functionals 



The phase behavior of a polydisperse fluid is determined by the chemical potential of 
the various species a, p = p(a,T,[p\), and the pressure, p = p(T,[p}), of the system. 
These functionals can be obtained from f(T, [p]) by the usual rules^l, suitably extended 
to a continuous mixture!, as: 

p(a,T,[p]) = 5f(T,[p])/5p(a) ; (1) 

p(T, [p})= J dap(a)p(a, T, [p]) - f(T, [p]) . (2) 

In what follows it will be convenient to split these quantities into an explicitly known ideal 
gaz part (id) and a generally unknown excess part (ex), for ex. f — fid+ fex, with: 

f id (T, [p]) = k B T J dap(a) {in (A 3 (<tM^)) - l} (3) 
p id (a,T, [p]) = k B T ln{AV)p(a)} (4) 

p id (T, [p]) = k B T J dap(a) (5) 

where k B denotes Boltzmann's constant, A(cx) the thermal de Broglie wavelength of species 
a, and the integrals over a extend over the whole domain of definition of p(cr). 



B. Phase coexistence conditions 



When for a given T, a parent phase of density Po{c), phase separates into n daughter 
phases of density pi(o) (i — 1, . . . ,n), the thermodynamic conditions of phase equilibriumS 
imply the equality of the pressures: 

p(T,[ Pl ])=p(T,[p 2 ]) = ...=p(T,[p n ]) (6) 

and of the chemical potentials: 
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Pi", T > [Pi]) = Vfa T > M = ■■■= P(&, T > \Pn\) (7) 

for each species a. It will be convenient to rewrite, Pj(cr) = Pjhj(a), where pj = J dapj(a) 
is the average density and hj(a) the normalized (viz. / dahj(a) = 1) polydispersity distri- 
bution of phase j (j = 0, 1, 2, . . . , n). The phase separation is moreover constrained by the 
conservation of the total number of particles of each species a: 

n 

^o(c) = ^Xihi^a) (8) 
i=\ 

where, Xi = Ni/N , is the ratio of the total number of particles (N) in phase i to the total 
number of particles in the parent phase (No). Moreover, the conservation of the total volume 
occupied by the parent phase can be written: 

n 

V = Y1 X i V i ( 9 ) 

where v j = 1/ Pj (j = 0, 1, 2, . . . , n). Finally, the normalization of the hj(a) in @ implies: 

1 = X> (10) 

i=i 

which expresses the conservation of the total number of particles. 

In principle, given T, p and ho(a), one has to solve the system of equations (|6|-[l0|) for 
the pi and hi(a) (i = 1, 2, . . . , n). Even when starting from a relatively simple f(T, [p]) this 
turns out to be a rather formidable task because (|]-0) are no longer algebraic equations (as 
would be the case for ordinary mixtures) but become here integral equations for the hi(a). 
In order to simplify the situation somewhat we will restrict ourselves here to two-phase 
coexistences only, i.e. to n = 2. In a polydisperse system n can in principle be arbitrarily 
large because the Gibbs phase rule! does not restrict the value of n in a system of infinitely 
many speciesi. In practice, however, a polydisperse system does rarely use this infinite 
number of thermodynamic degrees of freedom. Since these multiple-phase coexistences are 
expected to occur at low temperatures, the restriction to two-phase coexistences (n = 2) 
implies that the value of T should be chosen high enough (see below). Note that this is 
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consistent with our earlier restriction (see above) to fluid phases because for low temperatures 
some of the (multiple) fluid phases may have to compete with some of the solid phases. 

For a two-phase coexistence (n = 2) the number fractions X{ (i = 1,2) can be obtained 
from (|-[I0D as: 

_ _ V -Vl , . 

X\ — , X 2 — 

Vi -V 2 V 2 - V 1 

expressing the lever rule@, whereas h 2 (cr) can be eliminated by using moreover @: 

h 2 {a) = Aio(cr) H /ii((t) (12) 

where vq = l/po and h (a) are given parent phase data. To find V\ = 1/pi, v 2 = l/p2 and 
/ii(cr) we need three relations. To this end we use (f|) to rewrite (|7|) for n = 2 as: 

Pi(tf) = p 2 (ff)exp{ / 3A/i ex .(a,T, [pi,p 2 ])} (13) 

where /3 = l/k B T and : 

Ap ex (a, T, [pi, p 2 ]) = p ex ((J, T, [p 2 ]) - p e:c (a, T, [pi]) (14) 

with fi ex being the excess part of p(a, T, [p]). In terms of the ht(a), eq. ( |T3"| ) becomes: 

^i(o-) = h 2 (cr) Ai 2 (a,T; pi, p 2 ; [hi], [h 2 \) (15) 

where A\ 2 is a shorthand notation for A\ 2 = (vi/v 2 ) exp{/3Ap ex }. Eliminating /i 2 (cr) 
from (p~5|) by using (|12D one obtains finally: 

/ii((x) = h (cr) H (d, T; p , p 1; p 2 ; [/i ]) (16) 

where, H = {v 2 — vi) A 12 /{(v — V\) + {y 2 — v )Ai 2 }, and it is understood that h 2 (a) has 
also been eliminated from A 12 by using ([12]). Eq. ( |T^ ) is our first relation between p 1; p 2 
and /ii(a). Given T, po and ho(<j), eq. (p!6| ) can in principle be solved with respect to h\(a) 
for given pi and p 2 values. A first relation between pi and p 2 can then be found by solving 
moreover: 



1 = / dah (a) H(cr, T; p , Pi, P2, [h], [ho]) (17) 

which follows from (TI^| ) and the normalization of h\(a). Finally, the system of equations is 
closed by solving moreover (|J): 

p(T,[ Pl })=p(T,[p 2 ]) (18) 

expressing the equality of the pressures. It is obvious that solving a (two-phase) coexistence 
problem for a polydisperse fluid is much more involved than the corresponding problem for a 
monodisperse (single species) fluid. This is due mainly to the central integral equation fllBD 
which governs the change in the polydispersity or species distribution between the parent 
phase and the two coexisting phases (cf flT2|) ), a process usually called fractionation!. The 
extension of ( |i6| - |lTf ) for the case of a n-phase coexistence is straightforward. 

C. Thermodynamic stability conditions 

Of course, not for every prescribed {T,p (a)} will the corresponding phases exist or be 
thermodynamically stable. For this to be the case, the free-energy density, f(T, [p]), must 
remain a convex functional of p(a) = Po{<j), i.e. it must satisfy: 

f(T, [p + X5p}) < Xf(T, [p + 8p]) + (l- A)/(T, [p]) (19) 

for any, < A < 1, and for any change, Sp(cr) ^ 0, of the functional form of p{cr). If we 
consider only infinitesimal changes, Sp(a), eq. ( |T9| ) is equivalent to: 

<Etj(^-^)^/(t,[p]) (20) 

k=2 K - 

where 5 h f is the k-th functional variation of /: 

8 k f = J da,... J da k K k (a u ...a k ; T, [p]) 5p(a 1 )...8p{a k ) (21) 

and K k is an integral operator whose kernel is : 

K k (a 1 ,...a kl T,[p\) = —————— . (22) 

dp(ai)...dp{a k ) 
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Hence, the condition for stability with respect to infinitesimal changes (p(cr) — > p(a)+Sp(a)) 
reduces < 5 k f) to: 

0<5 2 /(T,[p]). (23) 

In a case of marginal stability there must exist at least one "fluctuation" 8p(o) 7^ 0, say 
5p (a), such that 5 2 f(T, [p]) = 0, or explicitly: 

J doi J da 2 K 2 (ai,a 2 ;T, [p]) 5p (a 1 )5p (a 2 ) = (24) 

for a given T and p(cr). Any solution Sp (a) of ( [2"4] ) will be called a critical fluctuation. 
For the system to remain stable with respect to these critical fluctuations 5p (cr) we must 
have according to (POfl, 5 3 f = (because the sign of S 3 f changes with the sign of Sp (a) 
while (Hp) does not fix the sign of Sp (a) ) and <5 4 / > 0: 

J dai J da 2 J da 3 K 3 (a±, cr 2 , cr 3 ; T, [p]) <5p ((Ti)5po(o'2)5po(o'3) = 0; (25) 

J dai J da 2 J da 3 J da 4 K 4 (at, a 2 , cr 3 , cr 4 ; T, [p]) 

.<5po(ai) 5p (cr 2 ) (5p (cT3) 5p (o- 4 ) > 0. (26) 



Note that in the present context the stability conditions (|23| - f26| ) imply stability with re- 
spect to changes in both the average density (Sp(cx) = 5ph(a)) and the composition 
(<5p(cr) = p5h(a)). The values of T and p(cr) for which ([24]) has a solution define a (general- 
ized or polydisperse) spinodal whereas those values for which fl24|) and ( p5| - |2~6|) are simultane- 
ously satisfied correspond to the critical states of the polydisperse fluid. For a monodisperse 
fluid these conditions reduce to well known results but for a polydisperse fluid they are much 
less well known and also much more difficult to study. For instance, for eq. ( |24T) to have a 
solution, the Fredholm determinani0 of the integral operator K 2 must vanish. This func- 
tional determinant however is equivalent to an infinite series of ordinary determinants^. As 
shown by Kincaid et al.0, in order to draw any conclusion from such a series of determinants 
one must have at its disposal a smallness parameter allowing one to truncate this series. For 



the weakly polydisperse fluids considered by Kincaid et al. 113 there is such a parameter but 
for the strongly polydisperse systems to be considered here there is none. Needless to say 
that in order to make progress we will need further simplifying assumptions. To this end we 
will introduce in the next section a model expression for f(T, [p]) which is simple enough to 
allow us to tackle the various technical problems encountered above in the study of phase 
equilibria in polydisperse fluids. 



III. THE POLYDISPERSE VAN DER WAALS FLUID 

An explicit expression for f(T, [p]), which is both relatively simple and physically sound, 
can be obtained from the van der Waals (vdW) theoryS. For a polydisperse fluid the vdW 
free-energy readsEH 

f(T, [p]) = k B Tj dap(a){H A3{ ^ a) ) - 1} 

+ lJdcrj da>V(a,a')p(a)p(a>) (27) 

where the hard repulsions between the particles are taken into account via the usual vdW 
excluded volume correct ionS, E[p]\ 

E[p] = 1 - J dav(a)p(a) (28) 

v(a) being the volume of a particle of species a, while the cohesion energy0 resulting from 
the interparticle attractions described by V^(r; a, a') is given in the vdW mean field approx- 
imation by the second term of (|27"D with: 

V(a,a') = J drV A (r;a,a'). (29) 

From (|27| ) and ([!]) we obtain: 

p(a,T, [p]) = ksT\n{ A3{ ^ a) } + fc B T^ J do' p{o>) 

+ J da' V(cr, a')p(a') (30) 
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while (H) yields: 



p(T, [p]) = ^J dop{a) + \jdojda> V(a, a')p{a)p{a'). (31) 

If D(a,a') denotes the contact distance between two particles of species a and a', we have 
Va(t; a, a') = when r < D(a, a') so that (^9|) can always be rewrittenSS 

Air 

V(a 1 a') = -e(a,a') Y (D(a 1 a')f (32) 

where e(cr, a') > characterizes the amplitude of the attractions. For simplicity we will use 
the Lorentz-Berthelot mixing rulesBi: 

D(a, a') = {D(a, a) + D(a', a')}/2 

e(a,a') = {e(a, ( r)e(a',a')} 1 / 2 (33) 

together with the additivity assumption, D(cr,a) = 2R(a), R(cr) being the radius of a par- 
ticle of species a, so that v(a) = (4n/3)R 3 (a). Choosing a = 1 as a reference species we 
may use -R(l) as length scale and e(l, 1) as energy scale. Our basic (dimensionless) poly- 
dispersity variable will hence be, R(a)/R(l), which will be denoted here, a = R(a)/R(l), 
so that a denotes here both a species and its associated polydispersity variable, R(a)/R{l). 
Below it will be convenient to consider various particular cases of (|32] - |33[) by introducing 
two indices, < / < 1 and < n < 1, and rewrite D(a,a) = a n D(l,l), v(a) = a 3n v(l), 
e(cr, a) = a 2l e(l,l) so that {/ = 0, n = 0} corresponds to the absence of polydispersity, 
{I — 1, n — 0} to polydispersity of the amplitude or interaction strength only, {/ = 0, 
n — 1} polydispersity in size only, while {I — 1, n — 1} corresponds to both amplitude and 
size polydispersity. Note that the amplitude polydispersityEl (/ ^ 0) includes the polydisper- 
sity which originates from the species dependence of the range of the interaction potential 
{Va{t; a, <j')). Finally, the present vdW-model is only a particular case of the general class of 
models considered by Gualtieri et al0. It has, however, the advantage that it can be easily 
related to a realistic interaction potential via Eqs. (E^) and (|32|). 
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IV. THE PARENT PHASE DISTRIBUTION 



To complete the description of our system we shall now specify, po(o) = p h (o), the 
distribution of the parent phases to be considered here. The polydispersity distribution 
ho(a) is constrained by the fact that it must be non-negative and normalized. We will 
assume moreover that the values of o = R(o)/R(l) are distributed continuously within 
< o < oo. Of course, very large values of o are unphysical but these will be given a very 
small weight by requiring that ho(o) decays with o in a manner which is sufficiently rapid 
for all the moments of h (o) } = J °° doo k h (o) , to exist. Moreover, we will restrict 
ourselves to monomodal distributions! so that the systems considered here are polydisperse 
generalizations of single component systems. Two often used candidatesiil are the Schulz- 
Zimm (SZ) distribution, }iq SZ \o) = co a exp(— bo), and the log-normal (LN) distribution, 
h,Q (o) = c exp{— a (\n(o/b)) 2 }. In each case the parameters {a, b, c } are determined 
by the normalization (m^ = 1), the mean value (mf ) and the polydispersity index, / = 
777,2 /( m i ) 2 > °f the distribution. In what follows we will take mf' = 1 so that the size 
of the reference species, R(l), is equal to the average value of R(o) in the parent phase 
distribution. Under these circumstances the SZ distribution can be written: 

hi SZ \o) = ^--o a - 1 -eM-^) (34) 

where T(a) is the Euler gamma function of argument < a < oo, a parameter which 
determines the inverse width of the SZ distribution, or equivalently, its polydispersity index, 
j(sz) _ -L _|_ \/ a _ moments of ([34]) are given by: 

(sz) _ J_ T(a + k) 
mk ~ a* r(a) { } 

which, when k is an integer, can be rewritten (for k > 2) as: 

-r ) = nti(i + -)- (36) 

a 

For the LN distribution we have similarly: 
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/TV)= 1 ^expU^l (37) 
(27rln/) 1/2 1 21n/ J v ' 

where j( iiV ) = 7 while the moments of (|37| ) are given by: 

m (LA0 = /fc (fc-l)/2 (38) 

In both cases we have hence mj, ' = mf^ = 1 and = I = 1 + 1/a. Note that in 
the monodisperse limit, {I — > 1 or a — > oo) both these distributions reduce to the Dirac 
distribution, <5(<r — 1), centered on the reference species (a = 1). On the contrary, when i" 
becomes very large (/ >> 1) the above distributions become very wide, increasing hereby 
the importance of the larger particles (see Fig.l). As stated above, the presence of very large 
particles is unphysical and a realistic situation should be described by distributions which 
strictly vanish for a larger than some maximum value. To avoid such problems, below we 
will consider only values of / in the range 1 < / < 2. For such values the results obtained 
from the full distributions will be similar to those obtained from the physically truncated 
distributions (with a slightly different /-value). 



V. THE SOLUTION METHOD 

In what follows it will be convenient to use dimensionless variables. Using e(l, 1) as 
energy scale and -R(l) as length scale, a = 1 being the reference species, we can introduce 
t = ksT I e(l,l), /I = /i/e(l,l) as reduced temperature and chemical potential while, if 
v (1) = ^-R 3 (l) represents the volume of the reference particle, for the reduced free-energy, 
pressure and average density we will use, respectively, / = fv(l)/e(l, 1), p = pv (l)/e(l, 1), 
1] = pv(l). The generalized moments (mfc(cx) = o k but k need not be an integer) of the 
polydispersity distribution h(a) will be written: 

rod 

m k[h] = / da mk(cr) h(o) (39) 

or, shortly, rrik = irik[h] for a general distribution h(a) and mj*^ = rrik[hj] when hj(a) 
represents the distribution of phase j (j = 0,1,2). Finally, as explained at the end of 
section [ill], V(cr, a') will be written: 
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( n n _|_ n-' n \ 3 

with 1/(1,1) = —8 e(l, l)f(l), where the value of / controls the amplitude polydispersity 
while the value of n controls the size polydispersity. In terms of these variables we can 



rewrite the excess part of eqs. (|30|), needed for (|14D, as: 

trj m m 3n (a) 



Vexiv, *> V, [h]) =-t ln(l - 77 m 3n ) + 



1-V m 3n 

'V {mi+zn{°)mi + 3mi +2n (a)m l+n + 3m l+n (a)m l+2n + mi(a)m l+3n } 



(41) 



while ( |3TD becomes: 



. . - tTJ TTZq > 

P\t, V, [h\) = 1 ?7 {mimi+3n + 3 m i+n m i+2n ) (42) 

1 - 7/ m 3n 



where it is seen that because of the polynomial character of (|40|) these expressions involve 
only a finite number of moments m^, {k = 0, 3n, Z, I +n, I + 2n, I + 3n}. Note that the exact 
number of moments depends on the values of I and n (for ex. when I — 0, k — {0, 1} or 
A; = {3n, / + 3n} represent one and the same moment, resp. k = and k = 3n). In the 
present context we can thus rewrite fllTf ) as: 

h x {a) = h (a) H(t, 770,771,772; {m k (cr)}, {m™}, {mf ] }) (43) 

where the {m k 2 } have been eliminated in favor of {m^. } and {m| } by using (|T^): 

(2) V 2 -V X ( o) V 2 -W (1) (AA s 

ml — ml, H ml.'. (44) 

Although eq. ([43]) can be solved directly, as was done in@, here we will take advantage of 
the structure of (4j) to transform this integral equation for h\{a) into a set of moment 



relations^: 

m k' )= / (icrm fc /(a)/7 (o-)^(t,^o,?7i,?72;{7r7fc(o-)},{m^ ) },{mi 0) }) (45) 



where both /c' and /c run through the finite set of moments {0, 3n, I, I + n, I + 2n, / + 3n}. 
Note that for fc' = 0, eq. (^) corresponds to (|I~7D while (|18D becomes here: 
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tyi m£ ] 2/ (i) (i) , o (i) (i) \ _ tr /2 mf ] 
- Tjr — T] x \m l m l+3n + 6 m l+n m l+2n ) — - t^t 

1 - Vi msn 1 - V2 m\n 

2/ (2) (2) . Q (2) (2) x 

-rj 2 {m\ X+3n + 3m^ n m^ 2n ) 



(46) 



together with the normalization condition ttiq = m^ 1 ' = rriQ = 1. For any given t, T]q 
and ho(cr) (from which the can be determined) the seven unknowns {?7i,^2,^3n ) m i > 
mi+ n i m i+2ni m i+3n} can t nen be obtained by solving the system of seven equations ( |45|) 
and fl4"6|) for {A;' = 0, 3n, I, I + n, I + 2n,l + 3n}. When this result is substituted in fl4"5|) 



we obtain h±(a) and from (O) we obtain then ^(c)- This then completely solves the two- 
phase coexistence problem for the present vdW-model. In the monodisperse limit we have 
mj^ = 1, for all k and j values, and reduce then to the usual vdW equations for 

the binodal of the reference species. To solve (f|5|-[|6|) we have used an iterative process by 
starting from an initial guess of the solution (usually the monodisperse case). When using 
moreover a globally convergent Newton- Raphson method0 a (convergent) solution of (|4l^46|) 



can easily be obtained even close to the critical points. The latter can hence be determined 
directly from the binodals but they can also be obtained by alternative procedures. One such 
procedure which is based on the stability analysis of section |II C| will be detailed in the next 
section whereas here we will focus on an alternative procedure based on the determination 
of the so-called cloud-point (C) and shadow (S) curves. These curves provide envelopes for 
the binodals. They can be obtained from ( |45}|46D by considering a situation of insipient 
phase separation whereby phase 1 is present only in infinitesimal amounts. Returning to 
section [11 B| this situation is seen to correspond to x\ — > or v% — > Vq with v\ and t>2 finite. 



From (|12|) and ([44]) it is seen that this implies h 2 (a) — > h (a) and mj^ — > m k °\ These curves 
are hence solutions of (cf. (|43| , |45|446|) ): 

h x {a) = h (a) A 12 {t,r] = 772,771; {m k (a)}, {m£>}, {m { £>}) (47) 



m£) = I dam k >(a)h (a) A 12 (t,r] = r] 2 ,r] l ]{m k (a)},{ml L, },{m^'}) (48) 
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trji mjp 



2/ (1) (1) , o (1) (1) \ 



2 ; (0) (0) .-, in) 1 1 1 1 

V ( m l m l+3n + 3 m l+n m l+2n) 



trjp m[, 0) 
1 (o) 

1 - V0 ™3n 
(0) ™(°) 



(49) 



where we took moreover into account in (|47|-48[) that, H — > A\ 2 , when v 2 — > vq (cf. (|T^- 
16|) . The solution to ( [4~7|]4"9l ) corresponding to the majority phase 2 yields then the so-called 
C-curve while the solution for the minority phase 1 yields the S-curve. 



VI. CRITICAL POINT ANALYSIS OF THE VDW FREE-ENERGY 

As explained in section [11 C| , a free-energy density f(T, [p]) will be globally stable if it 
remains a convex functional of p(a) for any change, p(a) — > p(a) + 5p(cr). It then must satisfy 
eq. (|l~9l). The same free-energy will be locally stable when (|2~0"D holds, i.e. when the first 
functional variation 5 k f(T, [p]) of f(T, [p]) which is non-zero corresponds to an even value 
of k > 2. We then have {S 2 f > 0} for a (ordinary) stable state; {5 2 f = 0, 5 3 f = 0, 5 A f > 0} 
for a (ordinary) critical state, {S 2 f = 5 3 f = 8 f = 5 5 f = 0, 8 6 f > 0} for a tricritical 
state, etc. Here we will limit ourselves to the equations determining the ordinary critical 



states of the vdW free-energy of section |T|. From (^7]) and (|22]) we obtain ((3 = 1/ksT] 
v(a n ) = v{a n )/E[p\): 

(3K 2 (a u a 2 ; T, [p]) = 5 ^ ~ ^ + p V (a u a 2 ) + [vfa) + v{a 2 )\ + p v{a x ) v(a 2 ) (50) 

ais t t r n 8(a x - a 2 )8(a 2 - a 3 ) _, , _, , . , 

/3A 3 ((Ti, a 2 , a 3 ; T, [p\) = — — — r h 2p v{ai) v{a 2 ) v(a 3 ) 

p{ai)p[a 2 ) 

+ [v(a 1 ) v(a 2 ) + v(a 2 ) v(a 3 ) +v(cr 3 ) v(ai)} (51) 

PA 4 ((Ti, CT 2 , <7 3 , CT 4 L , [p J = 2 r r ■; r 

P{°\) PM P{°3) 

+6p v(oi) v(cx 2 ) v(a 3 ) TJ(cr 4 ) + 2 [v(ai) v(a 2 ) v(a 3 ) + v(a 1 ) v(a 2 ) TJ(cr 4 ) 

+v(a 1 ) v(a 3 ) v(a 4 ) + v(cr 2 ) v(a 3 ) v(a 4 )] (52) 
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where it is seen that K4 is positive definite while K3 and K 2 can vanish because the volume 
v(a) of species a and E[p] of (|28| ) are positive while V (01,02) °f (2D * s negative. We 
now turn to eq. ( p^) which can be cast into an eigenvalue form0lil by rewriting Sp(a) as 
5p(cr) = p(a) ■ e(cr) where = ph(a) and e(<r) is as yet unknown. Using moreover the 
dimensionless variables introduced in section [V] we can recast (^|) and flSTf ) in the form: 

J daida 2 e(ai){(h(a 1 )h(a 2 )) 1/2 5(a x - a 2 ) + hfa) C 2 (a u a 2 \t,rj) h{a 2 )}e{a 2 ) = (53) 

where: 

C 2 (ai, a 2 ; t,rj) = A (m 3 „((Ti)mo((T2) + m (<7i)m 3n (cr 2 )) + A 2 m 3n (cri) m 3ri (cr 2 ) 

- S{mi +3n (o-i)m/(cr 2 ) + 3m i+2n (ai)mi +n ((j 2 ) 

+ 3m l+n (a l )m l+2n (a 2 ) + TO/(oi)m z+3n (c7 2 )} (54) 

together with the shorthand notations, A = rj/(l — r]m 3n ), and -B = 77/I For convenience 
we rewrite ([54]) briefly as: 



C 2 (cri,cr 2 ;t,rj) = ^ c kk >{t, rj)m k {ai)m kl {a 2 ) (55) 
fe.fc' 

where the c^/ can be deduced by identification of fl55j) with fl54[). From (^5j) it is obvious 
that C 2 (o"i, <T 2 ; t, r/) is completely embedded in the subspace spanned by the vdW-moments 
{m k (a); k = 0, 3n, Z, Z + n, Z + 2n, Z + 3n}. To solve ( |53"D it will hence be sufficient to restrict 
e(er) to the same subspace@, i.e. e(cr) = Z)fc e fc m A:( cr ) • A solution of the homogeneous 
equation (|53|) will hence exist provided the determinant of the matrix, M+M-C-M vanishes. 
We can rewrite this matrix as M_ -A, with A = L + M. and B_ = C ■ M. Here, I kk i = 5 kk , is the 
unit matrix, C kk > = c kk '(t,T]) the matrix of coefficents of (|55|), and M kk < = m k+k > a matrix 
involving the moments of h(a). (Note that m k (a)m k /(a) = m k+k i(a)). Since det\M kk /\ ^ 0, 
the equation of the spinodal (SP) can be written in various forms, e.g. : 

det\I kk i + B kk i \ = 0, B kk i = Yc kk »m k » +k f (56) 

k" 



whereas the corresponding solution of ([53]) reads (up to a proportionality factor): 
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e{a) = ^a k , k m k (a) (57) 



where a k < k is the co-factor (algebraic minor) of Ay k = 5j5l + B k r k (note that the r.h.s. of 
is in fact independent of the value of k'). Having found the spinodal (SP) we now turn to 
the stability condition (ST), eq. (p5[). Using fl5"T|) and writing Sp Q (a) = p(a)e(a), eq. fl25|) 
can be rewritten (E = E [p]): 

J da h(a) e 3 (a) = 3 (-|) 2 { J da h(a) e(a)}{ J da h(a) e(a) v(a)} 2 

+2(^) 3 {Jdah(a)e(a)v(a)} 3 (58) 

which on choosing k' = in (|57| ) and substituting the result into fl58| ) yields: 

J2 a oka ok >a ok ,<m k+k/+k// = 3A 2 (J2 k 

+ 2A 3 (J2 k a 0k m k+3n ) 3 (59) 



where, as before, A — rj/(l — rjm 3n ). Finally, in the vdW model eq. (|26| ) is always satisfied 
because K4 of ( p2|) is positive definite. 

In summary, those thermodynamic states {t,7], h(a)} for which 5 2 f > are stable while 
the marginally stable states (5 2 f = 0) obey the spinodal condition (p6[). The critical states 
are stable marginal states, i.e. they must obey both eq. (|^) and (|59|). Note that these 
equations are algebraic equations involving t, rj and moments of the type m k , m k+k i and 
m k+k i +k ii where {k, k', k"} cover the finite set of values {0, 3n, I, l + n, l + 2n, l + 3n} appearing 
in the excess free-energy. 

VII. SOME TYPICAL PHASE DIAGRAMS 

Having obtained the basic equations which govern the phase behavior of the present 
polydisperse vdW-fluid we now consider a few explicit examples. In order to completely 
characterize a thermodynamic situation we must specify (t, rjo), the parent phase distribution 
ho(a), and the values of (/, n) to be used in (^0[). In our study we have found no qualitative 
differences between the two types (SZ or LN) of parent phase distributions considered in 



section [IV]. All features found for one type can be found for the other type by changing 
the value of the polydispersity index /. On the contrary, the influence of the (l,n) values is 
much more pronounced, as we now explain. 



A. Interactions with size polydispersity only 

If we remove the amplitude polydispersity (I = 0) and keep only the size polydispersity 
(n — 1) of the interaction potential the present vdW model is governed by the following 
four moments {m , mi, m 2 , m 3 } , i.e. k = 0,1,2,3 (see section |V|). The equations (fl5|- 
f46|) governing the binodals constitute then a system of five equations which can be solved 
numerically and similarly for the C- and S-curves (see (|48|-49)), as explained in section |V]. 



Returning to (§^) it is seen that in the present case the effect of the polydispersity (I ^ 1) on 
the repulsions (or excluded volume) is stronger (rj — > rjm^) than its effect on the attractions 
(4rj 2 — > rj 2 (m 3 +3m 1 m 2 )) because the moments are increasing functions of / = 1 + 1/a > 1 
(cf. ( p6|) and fl38|)). This situation is hence unfavorable to phase separation when compared 
with the (strictly) monodisperse case (I = 1). This can be seen in the phase diagram 
shown in Fig.2. The result shown corresponds to a SZ parent phase of relatively modest 
polydispersity index I = 1.04 (a = 25). Note that at present^ it is current practice to 
consider colloidal dispersions with / < 1.05 as (approximately) monodisperse. It is seen in 
the figure that, when compared to the (strictly) monodisperse case (/ = 1), the coexistence 
region is shifted to lower temperatures. As indicated, there exists now a different binodal for 
each value of t]q, (instead of the unique binodal when 1=1). These binodals are truncated 
on the high-temperature side and fill the space between the cloud-point (C) and shadow (S) 
curves. The high-density phase is shifted towards lower densities whereas the low-density 
phase is only slightly shifted towards higher densities. Note also that the polydispersity 
distributions, hi(a) and ^2(0"), change with T along the binodals (see Fig. 3). The S-curve is 
situated in the interior of the C-curve and is tangent to the latter at the commun maximum 
of both curves. This commun maximum corresponds to a critical point where the two 
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phases become identical. The binodal through the critical point is the only untruncated 
binodal. The critical point itself is shifted to lower temperatures and lower densities when 
compared to the monodisperse case. When increasing the value of I all these shifts increase 
monotonically without modifying the overall aspect of the phase diagram. As far as we are 
aware, phase diagrams of this type have not yet been found experimentally but should be 
observable in suitably prepared colloidal dispersions. 



B. Interactions with amplitude polydispersity only 

If we remove the size polydispersity (n = 0) but keep the amplitude polydispersity (/ ^ 0, 
say / = 1) the present vdW model is governed by two moments {mo,mi}, i.e. k = 0, 1 and 
eqs. (|45| - |46|) reduce to a system of three equations. In this case the polydispersity does only 
affect the attractions (see, 4r) 2 — » rf ■ Am\, in eq. (|42|)), a situation which is favorable to 
phase separation. This is seen in Fig. 4 where the phase coexistence region is shown to be 
shifted to higher temperatures. Note also that the S-curve has moved now to higher densities 
and also moved partly outside of the interior of the C-curve. As a result the intersection 
of the C- and S-curves, i.e. the critical point, is encountered now below the temperature 
(t m ) corresponding to the commun maximum of the C- and S-curves. As seen in Fig. 4b this 
allows for a re-entrant behavior of the low density phase. The same behavior is found for 
any I > 0, whatever small. Phase diagrams of this type have been found for some polymeric 
systems^!. 



C. Size and amplitude polydispersity 

In a realistic situation we expect n = 1 and I ^ 0. This increases the number of moments 
(up to a maximum of six, five for 1 = 1) but it is easily verified that within the present 
vdW model the attractions are always more strongly affected by the polydispersity than 
the repulsions and will hence always favor the phase separation. Examples of the resulting 
phase diagrams are shown in Fig. 5. Note that the modifications introduced into the phase 
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diagram by the amplitude polydispersity (compare Figs. 2 and 5) could be used to probe the 
presence or absence of an amplitude polydispersity in the interaction between the colloidal 
particles being studied. 



VIII. POLYDISPERSITY INDUCED CRITICAL POINTS 

In the above we have limited ourselves to the case where only two phases coexist. We have 
seen that the phase diagram is considerably modified by the polydispersity, in particular in 
the region of the critical point. When the temperature is lowered, starting from the critical 
temperature, one expects to encounter a region where three, four, etc. phases coexist. 
Because only fluid phases are involved we expect the transition from the two-phase to the 
three-phase region to proceed through a second critical point. Since such a second critical 
point is absent from the ordinary monodisperse vdW fluid it can be termed polydispersity 
induced. In the present section we will describe how these polydispersity induced critical 
points can be found in a systematic way. Critical points can be obtained by following different 
routes, two of which have already been illustrated in section |V11| . It was seen there that 
critical points can be found by either looking for untruncated binodals or for intersections 
of the C- and S- curves. When the critical region is approximately known these routes are 
easily executed but this method becomes unpractical when no a priori knowledge about their 
location is available, as is the case here for all the polydispersity induced critical points. The 



method to be followed here will therefore be based on the stability criteria of section II C 



as adapted to the vdW case in section |VT|. When following this route the critical points 
can be found by looking for intersections of the spinodal (SP) and the stability (ST) curves. 
In the present context this amounts to look for thermodynamic states (t, 77) which, for a 
given h(a) = ho(a), satisfy both eq. (p6[) and fl59|). In Fig. 6 we show how the various routes 
coincide for the ordinary, high-temperature, vdW critical point. Some of the remaining 
polydispersity induced critical points obtained from the present route are shown in Fig. 7. 
It is seen there that when increasing the polydispersity index /, the number of polydis- 
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persity induced critical points gradually increases. It is also seen that at the same time the 
region of spinodal instability (5 2 f < for 5p(a) = Sp (a)) rapidly invades the high density 
region of the phase diagram. As a consequence, some of the polydispersity induced transi- 
tions and critical points could still be metastable with respect to the solid phases expected 
to become stable in this high density region. 

IX. CONCLUSIONS 

In the present study we have presented phase diagrams for fluids composed of spherical 
particles with a monomodal size distribution. These phase diagrams have been obtained 
on the basis of the van der Waals approximation for the free-energy of a polydisperse fluid. 
Interaction potentials with both size and amplitude polydispersity have been considered. 
It has been found that the largest modifications to the phase diagram of the polydisperse 
fluid, as compared to its monodisperse counterpart, result from the amplitude polydispersity 
of the interaction potential. These modifications should already be observable in colloidal 
dispersions with a relatively small polydispersity index / (e.g. 1.01 < I < 1.05). For larger, 
but still modest values of / (e.g. 1.1 < I < 1.5) a second, polydispersity induced, critical 
point is found signaling the onset of a three phase coexistence. When the value of / is 
increased still further, high-order coexistences are found but these are located in the high- 
density region of the phase diagram where the fluid phases considered here are expected to 
become metastable with respect to the solid phases. 
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Figure Captions 



FIG. 1. Two examples of the Schulz-Zimm (SZ) (solid lines) and log-normal (LN) (dotted 
lines) parent phase distributions, h Q (a), considered in this work (see fl34[ - j37|) ). The values 
of the polydispersity index (/) are, I = 1.04 and I = 4/3 (as indicated). It is seen that for 
the values of I considered here (1 < I < 2) these two distributions are rather similar. In 
particular, both are very small for the larger cx-values. Note also that the average value of 
a (mi = 1) does not coincide with the value of a for which h (a) reaches its maximum. 

FIG. 2. An example of a phase diagram (in (a) the reduced temperature (t = k B T/e(l, 1))- 
reduced density (r/ = pv{l)) plane and (b) the reduced pressure (p = pv(l)/e(l, l))-reduced 
temperature (t) plane) for the SZ parent phase distribution (/ = 1.04) shown in Fig.l and the 
vdW-model (I = 0, n = 1) with only a size polydispersity. For comparison the monodisperse 
phase diagram (/ = 1) (dashed line) has also been included. Shown are, the C-curve (full 
dots), the S-curve (open dots), and three binodals (full lines) for resp. 770 = 0.15, r] = 0.35 
and r/o = rjcrit = 0.2852 together with a tie line (dotted lines of (a)) corresponding to 
the upper temperature of coexistence (the lowest tie line shown belongs to the 770 = 0.15 
binodal). The critical point is indicated by a square (filled for / = 1.04 and open for / = 1). 
Note also that in the p — t diagram (b) the C- and S-curves coincide (with rji < 772 on the 
upper branch and 7/1 > 772 on the lower branch). 

FIG. 3. Evolution of the polydispersity distributions of the coexisting phases, hi(a) and 
/12(a), with the temperature along the critical binodal (770 = i] crit = 0.2852) of Fig.2. Shown 
are : hi(a) (full line, t = 0.9; dash-dot line, t = 1.12) and /12(c) (dotted line, t = 0.9; full 
dots, t = 1.12). 

FIG. 4. The same as Fig.2 but for the vdW-model (I — 1, n — 0) with amplitude polydis- 
persity only. The case shown corresponds to a SZ parent phase distribution with / = 1.02. 
The three binodals shown correspond to tjq = 0.25 (outer), 770 — 0.45 (inner) and the critical 
binodal 770 = r\ cr it = 0.346 (middle). Note that in the p — t diagram the critical point corre- 
sponds to the maximum of the pressure but not to the maximum of the temperature (t m ) 
for coexistence. The latter can also be seen in the t — rj diagram. This implies a re-entrant 
behaviour for t cr u <t <t m . 

FIG. 5. The same as Fig.2 but for the vdW-model (I — 1, n — 1) with both size and 
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amplitude polydispersity. The case shown corresponds to a SZ parent phase distribution 
with I = 1.01. The three binodals shown correspond to r] = 0.25 (outer), rj = 0.4 (inner) 
together with the critical binodal rj = r] crit = 0.3338 (middle). Note that this phase diagram 
is globally similar to the one of Fig. 4 but shifted to higher temperatures although the value of 
/ is smaller here. The phase diagrams are very sensitive to the total amount of polydispersity 
present, i.e. to the values of / and n. 

FIG. 6. The critical point shown in Fig. 5 as determined by three different routes : 1) the 
maximum of the untruncated critical (770 = i] crit ) binodal(B), 2) the intersection of the C- and 
S-curves and 3) the intersection of the spinodal (SP) and stability (ST) curves. Note that 
the SP-curve is tangent to the C-curve at the critical point (t crit = 1.2620, r] crit = 0.3338). 

FIG. 7. Evolution of the number of critical points for a vdW fluid with both size and 
amplitude polydispersity (I — 1, n — 1) and a LN parent phase distribution of increasing 
polydispersity ( (a): / = 1.1, (b): 7 = 4/3). Similar results can be obtained for the SZ 
distribution. The critical points (full dots) correspond to the intersection of the SP (full 
line) and the ST (dashed line) curves. Note that when I increases the number of critical 
points increases (one for 1 < I < 1.09, two for 1.09 < I < 1.5, etc.) while at the same 
time the spinodal region invades the high density portion of the reduced temperature (t)- 
average packing fraction (rj* = i]m 3 ) plane (here = J 3 is the third moment of the LN 
distribution). 
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Fig. 1 , Bellier-Castella et al., JCP 
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Fig. 2a , Bellier-Castella et al., JCP 




Fig. 2b , Bellier-Castella et al, JCP 




Fig. 3 , Bellier-Castella et al., JCP 




Fig. 4a , Bellier-Castella et al., JCP 




Fig. 4b , Bellier-Castella et al., JCP 




Fig. 5a , Bellier-Castella et al., JCP 
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Fig. 7 , Bellier-Castella et al., JCP 




